# Packages ----------------------------------------------------------------

library(haven)
library(tidyverse)
library(dplyr)
library(openxlsx)
library(zoo)
library(ggplot2)
library(tidyr)
library(flextable)
library(officer)
library(broom)
library(stargazer)
library(nnet)
library(Hmisc)
library(imputeTS)
library(marginaleffects)
library(mice)
library(miceadds)
library(sandwich)
library(clubSandwich)
library(modelsummary)
library(mitools) 
library(fixest)

# Load data --------------------------------------------------------------------

load("EFPC_data.Rdata")
full_data <- clean_data
rm(clean_data)

# Diagnostics -------------------------------------------------------------

## Missing data

missing_fraction <- colMeans(is.na(full_data))
sort(missing_fraction[missing_fraction > 0], decreasing = TRUE)
mean(missing_fraction)

# Imputation --------------------------------------------------------------------

full_impute <- mice(full_data, m = 20, method = "pmm", maxit = 10, seed = 12345)

# Split imputed data

full_impute.long <- complete(full_impute, "long", include = TRUE)
long_dem <- full_impute.long[which(full_impute.long$lex_binary == 1),]
long_aut <- full_impute.long[which(full_impute.long$lex_binary == 0),]
imp_dem <- as.mids(long_dem)
imp_aut <- as.mids(long_aut)

# Descriptive statistics --------------------------------------------------

## Full data ---------------------------------------------------------------

full_data |>
  ungroup() |>
  summarise(
    consult_1 = mean(consult == 1) * 100,
    no_consult_0 = mean(consult == 0) * 100
  )
## Regime type -------------------------------------------------------------

### Binary 

full_data |>
  group_by(lex_binary) |>
  summarise(
    consult_1_sum = sum(consult == 1), 
    no_consult_0_sum = sum(consult == 0),
    consult_1_percent = mean(consult == 1) * 100, 
    no_consult_0_percent = mean(consult == 0) * 100
  )

### Ordinal

full_data |>
  filter(!is.na(lexical_index_plus)) |>
  group_by(lexical_index_plus) |>
  summarise(
    consult_1_sum = sum(consult == 1), 
    no_consult_0_sum = sum(consult == 0),
    consult_1_percent = mean(consult == 1) * 100, 
    no_consult_0_percent = mean(consult == 0) * 100 
  )

## Consultation type -------------------------------------------------------

full_data |>
  filter(!is.na(lexical_index_plus)) |>
  group_by(lexical_index_plus) |>
  summarise(
    sum_pre = sum(stage_pre == 1), 
    sum_post = sum(stage_post == 1),
    sum_both = sum(stage_both == 1)
  )

## Region ------------------------------------------------------------------

full_data |>
  filter(!is.na(lex_binary)) |>
  group_by(region) |>
  summarise(
    consult_1_sum = sum(consult == 1), 
    no_consult_0_sum = sum(consult == 0),
    consult_1_percent = mean(consult == 1) * 100, 
    no_consult_0_percent = mean(consult == 0) * 100 
  )

## Decade ------------------------------------------------------------------

full_data |>
  mutate(decade = floor(year / 10) * 10) |>  
  group_by(decade) %>%                         
  summarise(consultation_rate = mean(consult))

# Logit models -------------------------------------------------------------

## All observations --------------------------------------------------------------

### Lexical Index of Democracy

full_logit_imp <- with(data = full_impute, 
                     exp = glm(consult ~ factor(lexical_index_plus) +
                                prior_const +
                                democratic_breakdown +
                                democratic_transition +
                                communist +
                                nationalist +
                                civilsoc +
                                protest +
                                civilwar +
                                avg_gdp_growth,
                                family = "binomial"
                                ))

### V-Dem polyarchy

full_logit_imp_alt <- with(data = full_impute, 
                            exp = glm(consult ~ polyarchy +
                                        prior_const +
                                        democratic_breakdown +
                                        democratic_transition +
                                        communist +
                                        nationalist +
                                        protest +
                                        civilwar +
                                        avg_gdp_growth,
                                      family = "binomial"))


## Democracies ---------------------------------------------------------------

dem_logit_imp <- with(data = imp_dem, 
                       exp = glm(consult ~ bodyctrl_perc +
                                   prior_const +
                                   democratic_transition +
                                   communist +
                                   nationalist +
                                   civilsoc +
                                   protest +
                                   avg_gdp_growth,
                                 family = "binomial"))

dem_logit_imp_interact <- with(data = imp_dem, 
                               exp = glm(consult ~ bodyctrl_perc*bodyctrl_new +
                                           prior_const +
                                           democratic_transition +
                                           communist +
                                           nationalist +
                                           civilsoc +
                                           protest +
                                           avg_gdp_growth,
                                         family = "binomial"))

## Autocracies ---------------------------------------------------------------

aut_logit_imp <- with(data = imp_aut, 
                      exp = glm(consult ~ performance + 
                                  partyinst +
                                  prior_const +
                                  democratic_breakdown +
                                  civilwar +
                                  communist +
                                  nationalist +
                                  protest +
                                  civilsoc +
                                  avg_gdp_growth,
                                family = "binomial"))


# Clustered SEs -------------------------------------------------------------

## SE function -------------------------------------------------------------

modelsummaryclus = function(model, list.imp){
  full_betas <- lapply(model, coef)
  full_vars <- lapply(model, FUN = function(x){vcovCL(x, cluster = list.imp[[1]]$country)}) 
  logit_se_results <- summary(pool_mi(full_betas, full_vars))
  ti <- data.frame(
    term = row.names(logit_se_results),
    estimate = logit_se_results[, "results"],
    p.value = logit_se_results[, "p"],
    std.error = logit_se_results[, "se"])
  mod <- list(
    tidy = ti,
    glance = NULL)
  class(mod) <- "modelsummary_list"
  return(mod)
}

## All observations -------------------------------------------------------------

full_impute_list <- miceadds::mids2datlist(full_impute)

full_logit_se <- with(data = full_impute_list, 
                            exp = glm(consult ~ factor(lexical_index_plus) +
                                        prior_const +
                                        democratic_breakdown +
                                        democratic_transition +
                                        communist +
                                        nationalist +
                                        protest +
                                        civilsoc +
                                        civilwar +
                                        avg_gdp_growth,
                                      family = "binomial"))

full_logit_results <- modelsummaryclus(full_logit_se, full_impute_list)

## Democracies -------------------------------------------------------------

### No interaction

dem_impute_list <- miceadds::mids2datlist(imp_dem)

dem_logit_se <- with(data = dem_impute_list, 
                      exp = glm(consult ~ bodyctrl_perc +
                                  prior_const +
                                  democratic_transition +
                                  communist +
                                  nationalist +
                                  civilsoc +
                                  protest +
                                  avg_gdp_growth,
                                family = "binomial"))

dem_logit_results <- modelsummaryclus(dem_logit_se, dem_impute_list)


### Interaction

dem_logit_interact_se <- with(data = dem_impute_list, 
                         exp = glm(consult ~ bodyctrl_perc*bodyctrl_new +
                                 prior_const +
                                 democratic_transition +
                                 communist +
                                 nationalist +
                                 civilsoc +
                                 protest +
                                 avg_gdp_growth,
                               family = "binomial"))

dem_logit_interact_results <- modelsummaryclus(dem_logit_interact_se, dem_impute_list)

## Autocracies -------------------------------------------------------------

aut_impute_list <- miceadds::mids2datlist(imp_aut)

aut_logit_se <- with(data = aut_impute_list, 
                      exp = glm(consult ~ performance + 
                                  partyinst +
                                  prior_const +
                                  democratic_breakdown +
                                  civilwar +
                                  communist +
                                  nationalist +
                                  protest +
                                  civilsoc +
                                  avg_gdp_growth,
                                family = "binomial"))

aut_logit_results <- modelsummaryclus(aut_logit_se, aut_impute_list)

# Figures -------------------------------------------------------------

## All observations -------------------------------------------------------------

full_plot_ord <- plot_predictions(
                full_logit_imp,
                conf_level = 0.95,
                type = "response",
                condition = "lexical_index_plus",
                vcov = ~country, 
                newdata = full_data) +
                xlab("Level of Democracy") +
                ylab("Predicted Probability") +
                scale_x_discrete(labels = c("Non-electoral \n autocracies", 
                                                     "One-party autocracies", 
                                                     "Multiparty autocracies \n (no elected exec.)", 
                                                     "Multiparty autocracies \n (elected exec.)",
                                                     "Electoral democracies", 
                                                     "Polyarchies")) +
                theme_bw() +
                scale_y_continuous(breaks = seq(0, 1, by = 0.1), 
                                   minor_breaks = seq(0, 1, by = 0.05)) +
                theme(axis.text.x = element_text(angle = 45, hjust = 1))

        
ggsave("full_plot_ord.tiff", 
       full_plot_ord, 
       width = 5, 
       height = 5, 
       dpi = 300, 
       compression = "lzw")

full_plot_cont <- plot_predictions(
                  full_logit_imp_alt, 
                  condition = c("polyarchy"),
                  type = "response",
                  conf_level = 0.95,
                  vcov = ~country, 
                  newdata = full_data) + 
                  xlab("Level of Democracy") +
                  ylab("Predicted Probability") +
                  scale_x_continuous(breaks = seq(0, 1, by = 0.1)) +
                  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, by = 0.1)) +
                  theme_bw()

ggsave("full_plot_cont.tiff", 
       full_plot_cont, 
       width = 5, 
       height = 5,
       dpi = 300, 
       compression = "lzw")

## Democracies -------------------------------------------------------------

dem_plot <- plot_predictions(
            dem_logit_imp, 
            condition = c("bodyctrl_perc"),
            type = "response",
            conf_level = 0.95,
            vcov = ~country, 
            newdata = dem_data) + 
            xlab("Seat Share of Largest Party/Coalition \n in Constitution-Approving Body") +
            ylab("Predicted Probability") +
            scale_x_continuous(breaks = seq(0, 1, by = 0.1)) +
            scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, by = 0.1)) +
            theme_bw()

ggsave("dem_plot.tiff", 
       dem_plot, 
       width = 5, 
       height = 5,
       dpi = 300, 
       compression = "lzw")

dem_plot_interact <- plot_predictions(
                    dem_logit_imp_interact, 
                    condition = c("bodyctrl_perc", "bodyctrl_new"),
                    type = "response",
                    conf_level = 0.95,
                    vcov = ~country, 
                    newdata = dem_data) + 
                    theme_bw() +
                    xlab("Seat Share of Largest Party/Coalition \n in Constitution-Approving Body") +
                    ylab("Predicted Probability") +
                    scale_x_continuous(breaks = seq(0, 1, by = 0.1)) +
                    scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, by = 0.1)) +
                    scale_color_discrete(name = "Newcomer status", 
                                         labels = c("0" = "Non-newcomer", "1" = "Newcomer")) +
                    scale_fill_discrete(name = "Newcomer status", 
                                        labels = c("0" = "Non-newcomer", "1" = "Newcomer"))

ggsave("dem_plot_interact.tiff", 
       dem_plot_interact, 
       width = 7, 
       height = 5,
       dpi = 300, 
       compression = "lzw")

## Autocracies -------------------------------------------------------------

aut_plot_partyinst <- plot_predictions(
                      aut_logit_imp, 
                      condition = c("partyinst"),
                      type = "response",
                      conf_level = 0.95,
                      vcov = ~country, 
                      newdata = aut_data) + 
                      xlab("Party Institutionalization") +
                      ylab("Predicted Probability") +
                      scale_x_continuous(breaks = seq(0, 1, by = 0.1)) +
                      scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, by = 0.1)) +
                      theme_bw()

ggsave("aut_plot_partyinst.tiff", 
       aut_plot_partyinst, 
       width = 5, 
       height = 5,
       dpi = 300, 
       compression = "lzw")


aut_plot_performance <- plot_predictions(
                        aut_logit_imp, 
                        condition = c("performance"),
                        type = "response",
                        conf_level = 0.95,
                        vcov = ~country, 
                        newdata = aut_data) + 
                        xlab("Performance Legitimation") +
                        ylab("Predicted Probability") +
                        scale_x_continuous(breaks = seq(-4, 2.5, by = 0.5)) +
                        scale_y_continuous(expand = c(0, 0), breaks = seq(0, 1, by = 0.1)) +
                        theme_bw()

ggsave("aut_plot_performance.tiff", 
       aut_plot_performance, 
       width = 5, 
       height = 5,
       dpi = 300, 
       compression = "lzw")

# Robustness checks ------------------------------------------------------------

## Linear probability model

### Full

full_lpm_se <- with(data = full_impute_list, 
                      exp = lm(consult ~ factor(lexical_index_plus) +
                                  prior_const +
                                  democratic_breakdown +
                                  democratic_transition +
                                  communist +
                                  nationalist +
                                  protest +
                                  civilsoc +
                                  civilwar +
                                  avg_gdp_growth))

full_lpm_results <- modelsummaryclus(full_lpm_se, full_impute_list)

### Democracies

### No interaction

dem_lpm_se <- with(data = dem_impute_list, 
                     exp = lm(consult ~ bodyctrl_perc +
                                 prior_const +
                                 democratic_transition +
                                 communist +
                                 nationalist +
                                 civilsoc +
                                 protest +
                                 avg_gdp_growth))

dem_lpm_results <- modelsummaryclus(dem_lpm_se, dem_impute_list)


### Interaction

dem_lpm_interact_se <- with(data = dem_impute_list, 
                              exp = lm(consult ~ bodyctrl_perc*bodyctrl_new +
                                          prior_const +
                                          democratic_transition +
                                          communist +
                                          nationalist +
                                          civilsoc +
                                          protest +
                                          avg_gdp_growth))

dem_lpm_interact_results <- modelsummaryclus(dem_lpm_interact_se, dem_impute_list)

### Autocracies

aut_lpm_se <- with(data = aut_impute_list, 
                     exp = lm(consult ~ performance + 
                                 partyinst +
                                 prior_const +
                                 democratic_breakdown +
                                 civilwar +
                                 communist +
                                 nationalist +
                                 protest +
                                 civilsoc +
                                 avg_gdp_growth))

aut_lpm_results <- modelsummaryclus(aut_lpm_se, aut_impute_list)
